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Abstract 


The  problem  of  propagating  acoustic  energy  in  the  ocean  subject  to  boundary  roughness  is  considered. 
Both  a  small  range-step  ” Monte  Carlo”  approach  in  which  the  acoustic  signal  is  scattered  from  an 
ocean  bottom  consisting  of  the  deterministic  bathymetry  with  a  stochastic  component  superimposed, 
and  a  scattering  kernel  method,  in  which  the  pressure  field  vector  is  analyzed  at  each  range-step  and 
modified  by  application  of  a  scattering  operator,  are  discussed.  Results  of  computations  using  the 
former  are  presented.  A  higher- order,  energy- conserving,  finite- element  parabolic  equation  model  is 
used. 
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1.1  Introduction 

The  importance  of  incorporating  interfacial  roughness  into  acoustic  propagation  models  is  clear.  This 
is  as  true  with  the  ocean  bottom  as  it  is  in  the  case  of  a  wind-roughened  sea  surface  or  the  under-ice 
canopy.  In  shallow  water  or  in  a  downward-refracting  environment  propagation  may  be  bottom- 
limited,  and  while  scattering  from  bottom  roughness  is  only  one  contribution  to  the  loss  of  energy 
from  the  propagating  signal,  it  is  known  to  be  important  in  some  situations.  The  same  is  true  in  the 
case  of  arctic  applications  or  sea  surface  roughness,  expecially  when  surface  interaction  is  strong. 
The  roughness  spectra  in  the  three  cases  are  clearly  distinct,  as  are  the  boundary  conditions, and  in 
the  case  of  the  bottom  there  are  many  complications  associated  with  transmission,  layering,  shear 
conversion,  etc.,  which  are  either  absent,  or  at  least  less  important,  in  the  other  cases. 

There  is  a  large  and  growing  literature  on  the  problem  of  scattering  from  finite  rough  surfaces. 
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growing  out  of  problems  in  both  the  electromagnetic  and  acoustic  domains,  and  especially  in  ocean 
acoustics  [1]-  [14].  Most  of  the  techniques  evolve  from  the  Helmholtz-KirchofF equation  and  involve 
different  ways  of  characterizing  the  field  on  the  boundary  surface.  Alternatively,  the  Lippmann- 
Schwinger  equation  has  been  employed,  and  occasionally  other  methods,  such  as  the  Green  s  function 
approach  of  DeSanto  [15] .  The  various  ways  of  attacking  the  problem  include  the  Rayleigh  hypothesis 
that  the  scattered  field  on  the  surface  can  be  represented  as  a  superposition  of  outgoing  plane  waves, 
even  though  that  is  known  to  be  valid  only  above  the  highest  point  on  the  surface.  The  extended 
boundary  condition  method  of  AVaterman  [2]  utilizes  the  fact  that  the  incident  and  scattered  waves 
must  cancel  in  the  region  below  the  interface  (in  the  case  of  an  impenetrable  boundary),  and  the 
extinction  method  forces  the  total  field  to  vanish  on  the  boundary  (in  the  Dirichlet  case).  Much 
of  the  effort  has  been  devoted  to  approximate  solutions  to  the  problem  of  rough  surface  scattering 
and  to  solving  the  Helmholtz  equation  ’’exactly,”  in  the  numerical  sense.  The  latter  includes  the 
work  of  Thorsos,  et  al[16],[17],  Fung  and  Chen  [19],  etc  [18],[20].  Other  work  has  dealt  with  the 
question  of  the  stability  and  uniqueness  of  the  solutions  to  the  Helmholtz-Kirchoff  integral  equation, 
especially  when  it  is  implemented  in  the  form  of  a  Fredholm  equation  of  the  first  kind  [21].  Further 
studies  have  invoked  perturbation  theory  or  the  Kirchoff  approximation  to  reduce  the  problem  to 
one  that  can  be  solved  explicitly.  The  extensive  numerical  work  has  drawn  a  rather  clear  picture  of 
these  approximations  and  their  ranges  of  validity.  Finally,  the  recognition  that  surface  roughness, 
especailly  that  of  the  bottom,  is  a  multi-scale  problem,  has  motivated  studies  of  composite-roughness 
models.  These  hybrid  or  two-scale  models,  which  divide  the  roughness  spectrum  into  a  large  and 
small-scale  part,  superimposing  a  small-scale  roughness  on  a  slowly  varying  surface,  suffer  from  the 
difficulty  of  determining  the  spatial  wavenumber  at  which  the  separation  between  large-  and  small- 
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scale  roughness  is  to  be  made,  and  seem  to  offer  more  promise  in  the  problem  of  scattering  from  the 
ocean  surface  than  in  the  case  of  the  ocean  bottom,  which  is  rough  at  "all”  scales,  i.e.,  the  relevant 
part  of  the  spatial  roughness  spectrum  spans  up  to  six  orders  of  magnitude. 

The  difficulty  with  all  of  this  work,  from  the  present  point  of  view,  is  that  it  sheds  little  light  on 
the  problem  of  incorporating  rough  interface  scattering  into  propagation  models,  and,  in  particular, 
PE  marching  models.  Indeed,  it  could  be  argued  that  the  divide  between  these  increasingly  elaborate 
attempts  to  exhaust  the  possibilities  of  describing  scattering  from  finite  rough  surfaces,  generally  in 
the  far  zone,  and  the  problem  of  roughness  in  propagation  modeling,  has,  if  anything,  grown,  rather 
than  diminished.  It  is,  in  any  case,  the  latter  problem  which  will  be  addressed  here,  although  it  will 
be  seen  that  in  an  approach  in  which  scattering  strengths  or  cross  sections  are  required  (see  below), 
they  would  be  obtained  from  the  work  summarized  above. 


1.2  Incorporating  Rough  Bottom  Scattering  into  a  Parabolic 
Equation  Model 

The  problem  of  interest  here  is  the  extension  of  these  ideas  to  the  propagation  of  sound  in  an 
acoustic  waveguide  with  rough  boundaries,  using  the  parabolic  equation  (PE).  In  fact  we  treat  only 
rough  bottom  scattering  because  of  the  difficulty  of  imposing  a  pressure-release  boundary  condition 
on  the  rough  (i.e.,  undulating)  ocean  surface,  but  in  principle  this  poses  no  essential  problem  for 
the  finite-element  PE  model,  and  Tappert  and  Nghiem-Phu  [22]  and  Thorsos,  Ballard,  and  Ewart 
[23], [24]  have  shown  how  to  deal  with  this  problem  in  a  split-step  model. 

The  simplest  approach  to  incorporating  some  of  the  effects  of  roughness  is  merely  to  lump  all 
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of  the  energy  which  is  scattered,  absorbed,  and  transmitted  at  an  interface,  and  treat  it  as  a  loss. 
In  a  more  sophisticated  implementation,  this  loss  will  be  angle-dependent,  as  in  the  treatment  of 
Moore-Head,  et  al  [25],  using  a  split-step  PE  model  in  which  the  angular-dependent  loss  is  applied 
in  wavenumber  space  via  an  FFT.  This  has  been  generalized  by  Dozier,  et  al  [26]  to  include  diffuse 
scattering  at  the  surface,  using  scattering  strengths  obtained  from  a  deterministic  model  of  the 
under-ice  canopy  in  terms  of  elliptical  bosses.  The  problem  with  this  approach  is  that  the  FFTs 
required  to  get  the  plane  wave  spectrum  at  the  surface  (or  bottom)  require  a  depth  window,  hence  a 
surface  layer,  with  the  result  that  the  scattering  takes  place  in  a  layer  rather  than  at  the  surface  and 
that  the  energy  emerges  from  the  surface  layer  down  range  from  where  it  encountered  the  surface. 
This  approach,  which  we  will  call  the  scattering  kernel  method,  is  also  being  investigated  in  detail, 
and  will  be  reported  upon  below.  At  the  moment,  we  briefly  disuss  its  implementation. 

Essentially  the  method  consists  of  decomposing  the  field  near  the  interface  into  a  plane  wave  spec¬ 
trum,  i.e.,  carrying  out  an  FFT  into  k-space,  then  operating  on  the  pressure  field  vector  in  k-space 
with  a  scattering  operator  matrix  elements  are  determined  by  calculation.  The  scattering  operator 
’’redirects”  the  energy  in  k-space  (angle  space)  at  each  range-step.  The  scattering  strengths  which 
go  into  the  scattering  kernel  may  be  determined  by  direct  numerical  integration  of  the  Helmholtz 
integral  equation  for  scattering  of  plane  waves  from  a  finite  surface,  by  application  of  the  Kircholf 
approximation  or  a  perturbation  method,  or  by  some  other  heuristic  method.  Since  the  scattering 
strength  depends  on  the  incident  and  final  grazing  angles,  the  roughness  parameters,  and  properties 
of  the  two  media,  the  storage  requirements  could  be  very  great.  A  better  solution  is  to  use  stored 
values  or  a  semi-empirical  expression  for  the  scattering  strength,  i.e.,  a  parameterization  of  the  scat¬ 
tering  cross-sections  so  that  the  process  is  reduced  to  a  look-up  table,  or  some  simple  closed-form 
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expression  e.g.,  Lambert’s  Law.  At  this  point  this  is  the  preferred  approach,  with  initial  work  us¬ 
ing  Lambert’s  Law,  although  we  will  ultimately  turn  to  an  expression  obtained  from  the  method  of 
small  perturbations  which  is  modified  to  give  approximately  the  correct  dependence  of  the  scattering 
strengths  on  frequency,  incident  and  scattered  grazing  angles,  and  roughness  parameters.  Assuming 
that  specular  reflection  is  already  correctly  treated  by  the  model,  it  must  be  extracted  from  the 
scattering  strength  before  it  is  applied  in  k-space. 

No  previous  attempt  has  been  made  to  apply  this  method  to  the  bottom,  where  the  interface 
conditions  are  in  general  more  complicated  than  at  the  ocean  surface.  Indeed,  the  method  makes 
a  realistic  treatment  of  the  ocean  bottom  (sediment,  basement)  difficult  or  impossible.  But  given 
the  uncertainty  principle,  i.e.,  the  need  for  a  finite  surface  layer,  there  seems  no  escape  from  these 
difficulties. 

One  is  entitled,  of  course,  to  ask  whether  the  scattered  energy  of  which  we  are  speaking  is 
important  at  al,  i.e.,  whether  it  is  important  to  propagate  it  down  range  (or  up  range  as  well,  as  in 
reverberation  studies),  rather  than  to  merely  treat  it  as  a  loss.  We  proceed  here  on  the  assumption 
that  this  is  indeed  the  case,  at  least  in  some  situations,  although  in  the  final  analysis,  it  will  be 
studies  such  as  this  which  will  answer  that  important  question. 

1.3  Small  Range-Step  or  ” Pull-Physics”  Approach 

The  only  remaining  realistic  approach  seems  to  be  to  directly  scatter  the  acoustic  energy  from  a  ran¬ 
domly  rough  interface  using  the  physics  built  into  the  model  itself.  Clearly  the  small-scale  roughness 
has  to  be  imposed  statistically,  and  the  essential  element  in  this  approach  is  the  small  range-step, 
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which  is  chosen  to  insure  that  the  model  ’’feels”  the  stochastic  roughness.  This  approach  we  call 
the  small  range-step  or  ’’full- physics”  approach.  Its  implementation  entails  generating  a  randomly 
rough  surface,  which  in  the  case  of  bottom  scattering  is  to  be  superimposed  on  the  deterministic 
bathymetry.  The  acoustic  energy  is  propagated  with  a  finite  element  PE  model,  subject  to  the  usual 
interface  conditions  on  the  boundary,  using  a  range-step  which  is  small  compared  to  a  wavelength. 
Through  the  interface  conditions,  the  model  will  reflect  energy  from  the  impedance  discontinuity 
and  the  results  become  very  like  application  of  the  Kirchoff  approximation,  in  which  the  scattering 
is  locally  specular.  Energy  conservation,  in  a  restricted  sense,  can  be  assured  by  requiring  continuity 
of  across  the  vertical  rises  of  the  stair-step  interface  [29].  The  question  of  energy  conservation 

has  been  addressed  by  Collins  and  Westwood  [29]  and,  more  generally,  by  Porter,  Jensen,  and  Ferla 

[30] . 

In  order  to  accurately  sample  the  statistical  roughness,  a  range-step  considerably  less  than  the 
acoustic  wavelength  is  used,  typically  TlOm  at  frequencies  of  25-125  Hz.  The  bathymetry  consists  of 
a  deterministic  part,  determined  ultimately  by  the  gridding  of  the  bathymetric  data,  with  a  Gaussian 
zero-mean  stochastic  component  superimposed  upon  it.  The  calculation  consists  in  marching  the 
finite-element  PE  solution  out  in  range  subject  to  the  impedance  matching  conditions  at  the  bottom, 
for  each  realization  of  the  statistically  rough  bottom.  The  final  transmission  loss  is  the  result 
of  averaging  over  the  N  realizations  of  the  bathymetry  (where  N  is  25-100),  each  with  the  same 
statistical  properties,  to  minimize  the  effect  of  statistical  fluctuations. 

Specifically,  if  we  decompose  the  complex  pressure  field  u  into  a  coherent  and  an  incoherent  part 

[31] 
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u  -<  u  >  +6u 


(1.1) 


where  clearly  <  Su  >=  0,  then  it  is  easily  shown  that 

<  |«|^  >=  I  <  «  >  |^+  <  |i5u|^  >  (1-2) 

where  the  terms  on  the  right  are  the  coherent  (specular)  and  incoherent  contributions  to  the 
received  intensity.  In  terms  of  transmission  loss, 


TLtoi=TLcoh+TLinc  (1-3) 

Clearly  TLtot  >  TLcoh- 

In  the  case  of  a  smooth  bottom,  some  energy  is  in  general  transmitted  across  the  interface,  while 
the  rest  is  specularly  reflected.  The  energy  transmitted  into  the  fluid  or  elastic  bottom  may  be 
refracted  back  into  the  water  column,  it  may  be  reflected  by  layered  structures  in  the  bottom,  etc.  If 
the  interface  is  rough,  however,  the  reflected  energy  will  be  scattered  into  directions  other  than  the 
specular,  in  general  smoothing  out  he  typical  pattern  of  convergence  zones  and  interference  effects. 


1.4  The  Parabolic  Equation  Model 

The  PE  model  used  is  a  finite  element,  higher  order,  energy-conserving  model  described  by  Collins 
and  Westwood  [29],  and  by  Collins  [32],  and  written  by  Michael  D.  Collins;  the  development  here 
follows  that  of  Ref.  32.  The  model  accurately  propagates  high  angle  energy  by  means  of  a  Fade’ 
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representation  of  the  square  root  operator  which  appears  in  the  parabolic  approximation  to  the 
Helmholtz  equation: 


+ -  0 


(1.4) 


For  two-dimensional  propagation,  the  (r,z)  equation  becomes,  in  the  far  field,  with  iheljy/r 
cylindrical  spreading  factor  removed, 


52$  a  1 5^  ,,2^  n 


(1.5) 


The  wavenumber  k  -  ujc  may  be  complex.  If  Eq.  (2.5)  is  factored  into  incoming  and  outgoing 
parts 


+  ikoK){-^  -  ikoK)^  =  0  (1.6) 

where  K  =  [I  +  k-^{p-§^^-§^  +  k^  -  kl)]^/'^. 

When  the  square  root  operator  in  the  outgoing  component  is  expanded  in  a  series  of  Fade’  approx- 
imants,  the  deceptively  simple  form  results: 


5$ 

—  =  ikoK^ 

or 


(1.7) 


l  +  a2i,M 

where  x  is  defined  by  /f  =  (1  +  The  Fade’  coefficients  are 


(1.8) 
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(1.9) 


2  .2 

=  2M  +  i*'"  2M  +  1 


==  """  2MTT 


(1.10) 


If  Eq.  (2.8)  is  written  (1  +  =  (ifcoaay-i.M^)^,  then  it  may  be  discretized  in  z  by  the 

Galerkin  method.  In  the  present  case,  the  radial  integration  is  carried  out  using  a  Crank-Nicholson 
algorithm. 

In  the  energy  conserving  PE  employed  here,  the  only  change  is  that  K  now  has  an  additional 
term  a  =  and  the  outgoing  equation  now  becomes 

^  =  ikoK^  (1-11) 

dr 


with  A"'  Fl  +  +  ik^  —  ib?!.  For  further  deails,  see  Collins  and  Westwood  [29]  or 

VYiljli  I  \  Ot  OZ  p  OZ  ®  J 

Collins  [32]. 


1.5  Bottom  Roughness  Spectra 

Many  studies[33],  [34]  have  shown  the  usefulness  of  characterizing  ocean  bottom  roughness  m  terms 
of  a  roughness  power  spectrum  which  exhibits  a  power  law  behavior.  Such  power  spectra  readily 
embody  the  fact  that  the  bottom  is  rough  over  five  to  six  orders  of  magnitude,  from  a  spatial  scale  of 
on  the  order  of  1  mm  to  1  km.  There  is,  in  fact,  a  natural  cut-off  at  low  spatial  wavenumbers,  where 
the  deterministic  part  of  the  bathymetry  takes  over,  and  at  very  high  spatial  wavenumbers  when  the 
roughness  scale  is  very  much  less  than  the  acoustic  wavelength  A.  This  is  not  the  place  to  explore 
the  largely  geophysical  question  of  the  power-law  parameters,  or,  indeed,  whether  a  simple  power 
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law  is  the  best,  or  even  an  adequate,  representation  of  the  statistical  properties  of  ocean  bottom 
roughness.  For  illustrative  purposes,  however,  we  will  discuss  the  effects  of  two  roughness  spectra, 
one  having  a  behavior,  with  a  rolboff  at  the  origin  [35],  the  other  being  a  Gaussian.  Specifically, 

we  use 


Wi{k)  =:  37r£crV[2(l  +  (1.12) 

W2{k)  =  a^£--;=zexp~^  ^  (1-13) 

v/tt 

Here  k  is  the  acoustic  wavenumber,  cr  the  rms  roughness,  and  £  the  correlation  length.  Although 
the  choice  of  roughness  power  spectrum  may  be  expected  to  be  important  in  certain  situations,  this 
is  not  the  case  in  the  examples  presented  here.  One  expects  that  the  frequency  dependence  of  the 
scattering  to  depend  on  the  roughness  power  spectrum  used,  since  that  determines  the  extent  to 
which  the  scale  of  the  roughness  is  commensurate  with  the  acoustic  wavelength. 


1.6  Application  to  Parabolic  Equation  Modeling:  Compu¬ 
tations 

In  this  work  we  employ  a  higher-order,  ’’energy-conserving”  PE  code  in  the  acoustic  environment 
shown  in  Figure  1;  the  bottom  parameters  are  given  in  Table  I.  The  bottom  roughness,  which  here 
consists  of  a  statistical  component  superimposed  upon  a  flat  bottom,  is  sampled  by  using  range-steps 
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of  1  m,  i.e.,A  /60  at  25  Hz  and  A  /20  at  75  Hz.  The  complex  pressure  fields  are  averaged  in  the 
manner  implied  by  Eqs.  (2.1)-(2.3),  using  N  realizations  of  the  statistical  roughness,  to  obtain  the 
total  and  coherent  intensities,  and  the  corresponding  transmission  losses.  In  general,  the  number 
of  realizations  is  made  as  large  as  is  computationally  feasible,  e.g.,  N=100,  since  the  statistical 
fluctuations  go  like  ^  ,  but  ordinarily  one  can  get  meaningful  results  for  N>  25.  The  realizations 
are  generated  by  a  Fourier  method  described  by  Thorsos.[17] 

1.7  Results  of  Numerical  Studies 

The  implementation  of  the  program  described  above  is  quite  straightforward,  though  computation¬ 
ally  intensive,  and  is  being  carried  out  on  the  CRAY  YM-P  at  Stennis  Space  Center,  MS.  Most  of 
the  results  presented  are  for  the  k~^  power  spectrum,  for  which  a  typical  realization  of  the  bottom, 
with  rms  roughness  cr  of  20m  and  correlation  length  t  of  500m,  is  shown  in  Figure  2.  In  the  test 
case  described  here,  the  acoustic  field  for  a  smooth  bottom  {a  =  0)  bottom  is  as  shown  in  Figure  3. 
The  75  Hz  source  is  at  150m  and  the  bottom  is  at  3000m;  again,  refer  to  Figure  1  and  Table  I  for 
the  acoustic  environment. 

The  effects  of  bottom  roughness  can  be  immediately  seen  in  Figures  4  and  5  in  which  an  rms 
roughness  of  cr  =  20m  and  correlation  length  I  =  500m  have  been  used  (this  value  of  t  is  used 
in  all  the  figures).  In  particular,  comparision  of  the  total  with  the  coherent  intensity  shows  that 
much  of  the  reflected  energy  is  scattered  diffusely;  the  incoherent  intensity  (Figure  6)  shows  this  fact 
dramatically.  Taking  the  receiver  to  be  at  150m,  one  sees  that,  especially  between  10  and  20  km, 
the  coherent  transmission  loss  is  lower  by  as  much  as  15  dB  (e.g., compare  TLcoh  in  Figure  7  with 
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(T  =  0  in  Figure  8)  and  that  beyond  the  point  where  the  sound  first  interacts  with  the  bottom  (10 
km)  [36],  TLcoh  is  always  considerably  less  than  TL(o«. 


SOUND  SPEED  (M/SEC) 


Figure  1.  Sound  speed  profile  in  the  water  column.  The  sediment  sound  speed  is  1700  m/sec. 
See  also  Table  1. 
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Table  1.1:  Acoustic  Parameters.  The  frequency  is  25  or  75  Hz 


Co  =  1500msec  ‘ 

=  0.0 

Q:jed(3000)  =  0.5 
a,ed(5000)  =  10.0 

npade'  —  2 
Zs  =  150m 

The  total  transmission  loss  in  Figure  8  exceeds  the  coherent  transmission  loss  by  as  much  as  20  dB 
in  places,  especially  beyond  10  km,  where  energy  from  the  first  bottom  bounce  reaches  the  surface. 
The  difference  then  settles  down  to  5-10  dB. 


Figure  2.  Example  of  a  single  realization  of  the  bottom  roughness,  with  <r  =  20m  and  £  =  500m. 
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DEPTH  (M) 


When  the  total  transmission  loss  for  <r  =  0  (smooth  bottom)  and  cr  =  20m  are  compared  (Figure 
8),  it  is  seen  that  beyond  10  km  the  transmission  loss  curve  is  somewhat  smoothed-out  by  the 
roughness  and  intensity  levels  are  generally  lower  by  a  few  dB.  Of  course  the  effects  are  smaller  for 
smaller  rms  roughness  (Figure  5). 


Figure  3.  The  acoustic  intensity  for  a  smooth  bottom  for  the  acoustic  environment  of  Figure  1 
and  Table  1.  The  frequency  is  75  Hz.  Intensities  are  plotted  from  a  TL  of  40  dB  to  lOOdB,  from 
white  to  black. 
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RANGE  (M) 


Figure  4.  Total  intensity  <  |up  >  for  rms  roughness  cr  =  20r7i  and  correlation  length  t  —  500m. 


RANGE  (M) 


Figure  5.  Coherent  intensity  |  <  u  >  ^  for  the  same  conditions  as  Figure  4. 
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RANGE  (M) 


Figure  6.  Incoherent  intensity  <  > 


Rongo  (m) 


Figure  7.  Total  transmission  loss  TLtot  and  Coherent  transmission  loss  TLcoh  ^t  75  Hz  for  source 
and  receiver  at  150m.  In  this  figure,  cr  =  20m  and  I  =  500m. 
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Ron9«  (m) 

Figure  8.  Total  transmission  loss  TL<„,  at  75  Hz  for  ct  =  0  and  a  =  20m  for  source  and  receiver  at 
150m. 

1.8  Frequency  Dependence 

Extensive  calculations  over  the  frequency  range  of  25-125  Mz  reveal  a  frequency  dependence  to  the 
scattering,  as  might  be  expected.  Although  this  is  difBcult  to  illustrate,  since  the  difference  in  the 
acoustic  fields  at  different  frequencies  is  not  by  any  means  due  to  scattering  alone,  it  appears  that  the 
effect  of  an  rms  roughness  a  =  20m  is  measurably  greater  at  100  Hz,  where  the  acoustic  wavelength 
is  on  the  order  of  the  rms  roughness,  than  at  25  Hz,  where  it  is  much  larger.  We  note  in  Figure 
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9  the  difference  between  the  total  transmission  loss  for  the  Gaussian  and  k  ^  power  spectra.  The 
differences  are  greatest  at  low  grazing  angles,  i,e.,  long  range,  as  would  be  expected. 

1.9  Discussion 

The  results  obtained  to  this  point  suggest  that  the  small  range-step  ’Tull-physics”  approach  gives  a 
plausible  description  of  the  specular  and  diffuse  scattering  from  a  rough  ocean  bottom.  Attention  is 
now  being  given  to  the  important  step  of  providing  quantitative  validation  of  these  results,  both  in 
a  real  acoustic  environment,  and  theoretically.  In  the  the  latter  case,  one  can  attempt  to  construct 
an  artificial  isovelocity  environment  in  which  the  bottom  is  smooth  except  for  a  scattering  patch 
at  some  range  R,  so  that  approximate  comparisions  can  be  made  between  the  PE  propagated  field 
and  numerical  simulations  using  the  method  of  moments.  It  is  also  useful  to  test  the  results  against 
the  results  of  a  normal-mode  propagation  code,  e.g.,  KRAKEN-C  [37],  in  which  roughness  is  taken 
into  account.  Finally,  one  can  attempt  intermodel  comparisons  between  the  two  general  methods  of 
attack  described  here.  Extension  of  the  small  range-step  approach  is  being  made  to  three-dimensional 
propagation,  where  a  2-D  roughness  power  spectrum  is  used. 
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Figure  9.  Dependence  of  the  total  transmission  loss  {TLtot)  on  the  roughness  power  spectrum  W(k) 
for  the  same  conditions  as  Figure  7.  See  Sec.  5  for  details. 


2.1  Scattering  Kernel  Method 

This  approach  offers  an  interesting  alternative  to  the  ’’brute-force”  small  range-step  metheod  which, 
in  principle,  does  not  entail  a  significant  increase  in  computational  overhead  compared  to  a  conven¬ 
tional  PE  calculation.  The  scattering  kernel  approach  consists  in  decomposing  the  pressure  field  near 
the  interface  (here,  the  bottom,  but  in  princple  the  surface  as  well)  in  vertical  wavenumber  space 
and  operating  on  the  field  with  a  scattering  operator  before  reconstructing  the  field  in  coordinate 
space.  Because  it  depends  on  obtaining  the  plane  wave  spectrum  at  the  interface-in  practice,  in  a 
layer  adjacent  to  or  spanning  the  boundary-two  FFTs  per  range  step  are  required.  The  need  to  keep 
the  surface  layer  thin  so  that  the  method  can  be  applied  in  shallow  water  and  the  desire  to  keep  low 
angle  rays  from  travelling  long  distances  in  the  surface  layer  requires  that  the  number  of  points  in 
the  EFT  window  be  severely  limited  (e.g.,  16,  32  points),  hence  the  angular  resolution  is  degraded. 
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If  the  scattering  kernel  is  computed,  in  the  near  zone,  at  each  range  step,  from  the  current  geometry 
and  bottom  properties,  then  the  approach  is  indeed  computationally  very  intensive.  If,  on  the  other 
hand,  some  simple  parameterization  of  the  scattering  from  the  bottom  is  used-the  present  work  has 
thus  far  employed  (for  simplicity)  only  a  Lambert’s  law  representation  of  the  scattering-then  the 
problem  is  much  less  severe.  Indeed,  all  of  the  manipulations  required  to  obtain  the  scattering  ma¬ 
trix,  which  depending  on  the  approach,  may  involve  diagonalizing  an  exponentiated  matrix  [26],  can 
be  done  once,  before  the  solution  is  marched  out  in  range,  or  at  worst,  only  when  the  deterministic 
bathymetry  changes. 

This  technique  was  first  proposed  by  Moore-Head,  Jobst,  and  Holmes  [25]  in  a  split-step  PE 
model,  in  order  to  introduce  angle-dependent  surface  loss.  The  generalization  by  Dozier,  Hanna,  and 
Pearson  [26]  introduced  a  scattering  kernel  which  redistributed  the  energy  incident  upon  the  surface, 
in  k-space.  Further  work  on  this  approach  has  been  reported  by  Schneider  [27].  All  of  the  physics 
is  in  the  construction  of  the  scattering  operator,  which,  however,  can  be  approximated  in  various 
ways  to  save  computational  resources.  These  include  adopting,  in  modified  form,  appropriately 
parameterized  and  scaled,  the  closed-form  expressions  for  scattering  in  the  perturbation  method,  or 
storing,  in  a  look-up  table,  the  results  of  numerical  scattering  experiments.  This  is,  in  any  case,  a 
matter  of  detail. 

The  method  proceeds  as  follows  [25]-  [27]  :  Given  the  complex  pressure  field  u(r,z),  multiplied  by 
a  smoothing  factor  w(z)  which  helps  prevent  aliasing  in  the  FFT  (typically  cos^  form  which  tapers 
to  0  at  the  upper  edge  of  the  enlarged  depth  window  of  thickness  2d,  with  R.,  the  cycle  distance  in 
the  layer  given  by  being  the  grazing  angle),  a  discrete  Fourier  transform  is  carried  out  into 

vertical  wavenumber  space.  The  field  vector  in  k-space  is  then  operated  upon  by  a  scattering  matrix 
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which  transforms  the  field  according  to  some  model  of  the  scattering  process,  and  finally  an  inverse 
FFT  restores  the  pressure  field.  The  procedure  can  be  represented  as  follows: 


rP{r,k)  =  FFT[w{z)U{r,z)] 

From  which  the  scattered  field  <i>{r,  k)  in  angle-space  is  found: 


(2.1) 


4>{r,k\)  =  ^^S{ki,k2)'4>{r,k2),  (2-2) 

where  S  is  the  scattering  operator  which  may  also  include  bottom  loss,  etc.  Finally,  the  pressure 
field  u  I  (r,z)  is  recovered: 


u\r,z)  =  FFT-^[ip(r,k2)]  (2-3) 

Dozier,  et  al  [26]  worked  with  the  intensity  |ti|2  rather  than  the  pressure,  and  showed  that, 
properly  formulated,  the  theory  guarantees  that  the  energy  will  not  grow  as  the  loss  and  scattering 
functions  are  applied.  In  this  case  the  critical  quantity  is  exp{SAr)  where  Ar  is  the  range-step  and 
S  is  the  scattering  kernel.  Although  S  is  not  precisely  the  same  operator  as  used  above,  it  is  again 
obtained  from  a  model  of  the  scattering  process  (in  Dozier,  et  al,  from  a  deterministic  model  of  the 
under-ice  surface  based  on  cylindrical  bosses);  in  either  case  any  loss  terms,  which  appear  in  the 
diagonal  elements  of  S,  represent  loss  per  range  step  over  the  cycle  distance  R.  In  the  formulation 
of  Dozier,  et  al,  the  phase  of  the  field  is  recovered  from  that  of  the  original  complex  pressure  field 
u(r,z). 
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There  are  some  good  objections,  of  course,  to  this  approach.  That  one  is  on  shaky  ground, 
numerically,  is  obvious.  By  modifying  the  PE  field  at  each  range-step,  one  is  in  danger  of  violating 
the  assumptions  under  which  the  parabolic  equation  was  derived.  More  serious  is  the  fact  that  the 
problem  is  no  longer  well-posed,  and  numerical  instability  is  the  likely  result.  At  the  same  time, 
it  does  seem  to  be  possible  to  carry  out  this  program  in  practice,  with  great  care,  as  the  work  of 
Moore-Head,  et  al,  Dozier,  et  al,  Schneider,  and  the  author  [28]  seems  to  show.  Nonetheless,  it  is 
not  clear  that  a  model  which  introduces  an  artificial  bottom  layer  when  the  real  bottom  has  its  own 
complicated  structure,  shear  conversion,  etc.,  is  desirable.  For  this  reason  the  method  is  likely  to 
have  more  application  to  the  ocean  surface. 

Thus  the  usefulness  of  this  approach  remains  to  be  established.  It  does  offer  the  possibility  of 
incorporating  bottom  roughness  into  something  other  than  a  research  code.  We  have  implemented 
the  method  in  a  finite-element  PE  code  using  a  scattering  kernel  obtained  from  a  simple  Lambert’s 
law  description  of  the  scattering.  Not  surprisingly  the  method  has  proved  to  be  susceptible  to 
aliasing  and  numerical  instability  so  that  it  is  not  always  easy  to  determine  which  results  are  real 
and  which  are  numerical  artifacts. 

The  spectral  decomposition  method  of  Gilbert  and  Evans  [38]  offers  an  alternative  war  of  finding 
the  plane  wave  spectrum  at  the  interface,  but  has  its  own  drawbacks,  and  in  the  end  all  methods 
for  finding  the  PWS  ultimately  come  up  against  the  uncertainty  principle.  As  mentioned  earlier, 
the  method  is  not  well-suited  for  shallow  water,  because  of  the  artificial  bottom  layer,  or  for  low 
frequencies,  where  the  layer  might  be  less  than  one  wavelength  thick.  It  is  certainly  not  well-suited 
to  treating  the  bottom  realistically. 
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3.1  Three-Dimensional  Parabolic  Equation  Modeling 

Given  the  success  of  the  ensemble- averaged  parabolic  equation  (PE)  treatment  of  scattering  from 
rough  surfaces,  it  is  of  interest  to  treat  three-dimensional  acoustic  propagation  in  the  ocean  waveguide 
subject  to  boundary  roughness  using  the  same  method.  The  presentd  study  focused  on  the  ocean 
bottom,  for  which  the  roughness  was  described  using  a  2-D  roughness  power  spectrum.  The  final 
acoustic  field  is  be  the  result  of  averaging  a  large  number  of  PE  runs  with  different  realizations  of 
the  statistically  rough  surface.  The  problem  is  computationally  intensive  but  naturally  amenable  to 
parallelization. 

One  of  the  central  problems  in  ocean  acoustic  propagation  is  that  of  incorporating  interfacial 
roughness  into  current  propagation  models  in  a  realistic  way.  At  the  air-water  interface  this  roughness 
will  consist  of  the  wind-roughened  sea  surface  or  the  under-ice  canopy.  In  the  case  of  the  seafloor,  the 
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roughness  includes  not  only  the  water-sediment  interface,  and  the  basement,  but  bottom  sediment 
layering  as  well.  Heretofore,  most  of  the  work  which  has  incorporated  boundary  roughness  has  done 
so  by  means  of  a  heuristic  reflection  coefficient  and  few  applications  have  been  to  parabolic  equation 
models.  Because  of  the  somewhat  greater  complexity  of  imposing  the  pressure-  release  boundary 
condition  on  the  rough  ocean  surface  (absent  ice),  we  shall  concentrate,  in  this  study,  on  the  ocean 
bottom,  and  in  particular  the  water-  sediment  interface. 

Of  the  two  approaches  described  above,  both  previously  explored  by  the  principal  investigator 
[39]-[41],  and  employed  in  Parts  1  and  2  of  this  work,  namely  the  ensemble-averaged,  small  range- 
step  and  the  scattering  kernel  methods,  it  is  the  former  that  was  applied  to  the  problem  of  3D 
acoustic  propagation,  using  a  two-dimensional  roughness  spectrum.  This  extension  is  motivated  by 
the  success  of  the  previous  studies  and  the  increasing  interest  in  true  three-dimensional  propagation 
in  ocean  acoustics.  More  generally,  the  importance  of  rough-bottom  scattering  lies  in  the  need  to 
accurately  predict  foward  propagating  acoustic  signal  levels  using  a  PE  model,  and  in  the  importance 
of  active  acoustic  applications  in  which  the  PE  field  is  propagated  backward  to  the  source.  Although 
it  is  adequate  in  most  situations  to  make  the  approximation  of  2-D  propagation,  that  is  not  true 
when  the  properties  of  the  water  column,  notably  the  sound  speed,  and  the  character  of  the  ocean 
surface  or  bottom  vary  rapidly  with  azimuth. 

Although  there  are  no  problems  in  carrying  out  this  program  in  principle,  there  are  two  major 
issues,  aside  from  those  which  arise  from  the  physics  of  the  problem.  The  first  is  the  need  to  optimize 
existing  PE  codes  for  the  present  purpose  or  to  develop  an  original  code  which  takes  advantage  of  the 
opportunities  to  use  parallel  processing.  The  small  range-step  approach,  which  ensemble-averages 
the  results  of  a  large  number  of  PE  runs  using  different  statistically  similar  realizations  of  the  rough 
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interface,  is  naturally  amenable  to  parallelization.  The  second  issue  has  to  do  with  the  2D  roughness 
spectrum.  Although  it  is  not  a  mojor  goal  of  the  proposed  study  to  address  the  issue  of  the  best 
representation  of  real  bottom  roughness  power  spectra,  the  effect  of  different  representations  of  the 
statistical  properties  of  the  bottom  has  been  studied,  as  well  as  the  relation  between  power  spectrum 
and  acoustic  frequency  dependence. 

This  final  stages  of  this  work  placed  great  demands  upon  the  computational  resources  available 
to  the  principal  investigator.  The  ensemble-averaged  small  range-step  method  is  inherently  compu¬ 
tationally  intensive,  due  both  to  the  nature  of  the  computational  grid,  the  need  to  average  a  large 
number  of  PE  runs  to  reduce  the  effects  of  statistical  fluctuation.  Its  use  in  3-D  PE  computations 
means  execution  times  at  least  an  order  of  magnitude  greater  than  in  the  2-D  case.  This  turns  out 
to  be  manageable,  however,  since  in  the  latter  case  the  acoustic  field  can  be  propagated  out  to  a 
range  of  50  km  using  aim  range-step  and  25-f  realizations  of  the  bottom  using  30  minutes  of  CRAY 
Y-MP8  time  without  parallelization. 

3.2  Three-dimensional  Parabolic  Equation 

This  section  describes  the  initial  phases  of  the  generalization  of  the  work  on  scattering  from  rough 
interfaces  in  two  dimensions,  using  a  roughness  power  spectrum  which  is  one-dimensional,  to  the 
case  of  full  3-D  propagation.  We  proceed  as  follows,  following  Collins  and  Chin-Bing  [42]; 

Starting  from  the  Helmholtz  equation  in  three  dimensions. 
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(3.1) 


K^P  =  0 

^  r  dr  dO'^  p  dz  dz 

Factoring  into  incoming  and  outgoing  parts  leads  to  an  outgoing  solution  of  the  form 

dP 


dr 


=  ikoQP 


(3.2) 


where 


5  =  +  +  (3,3) 

Here  K  is  the  complex  wavenumber.  If  the  rapidly  varying  plane  wave  e'*'*'"*  is  factored  out  of  P, 
the  approximations  of  Ref.  (4)  lead  to 


dP 
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(3.4) 


(3.5) 


(3.6) 


r^kl  dd^ 

The  result,  Eq.  (15)  is  a  wide-angle  3-D  parabolic  equation.  In  Ref.  (42),  this  equation  is  solved 
using  the  method  of  alternating  directions,  Galerkin  discretization,  Crank-Nicholson  integration  in 
r,  and  centered  differences  in  6.  The  approach  of  Lee,  et  al  [43]  to  the  solution  of  the  equivalent  of 
Eq.  (15)  is  to  expand  Q  somewhat  differently  and  to  use  different  numerical  techniques. 
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3.3  Numerical  Implementation 

The  most  direct  route  to  implementation  of  a  3-D  PE  scattering/  propagation  model  proved  to  be 
the  use  of  an  existing  three-dimensional  PE  code,  written  by  Michael  Collins  [42],  This  finite-element 
3-D  code  was  used  to  generate  the  three-dimensional  complex  pressure  fields  which  are  needed  to 
carry  out  the  ensemble  averages  described  above,  using  a  statistically  rough  seafloor  (in  general 
superimposed  on  whatever  the  deterministic  bathymetry  may  be)  obtained  from  a  2-D  roughness 
power  spectrum.  In  spite  of  this  fact,  considerable  modification  of  the  code  was  necessary,  and  it 
had  to  be  made  part  of  an  emsemble-averaging,  small  range-step  code. 

Although  development  of  the  scattering/propagation  code  is  straightforward,  the  demands  on 
computer  resources  are  naturally  quite  large.  In  the  two-dimensional  case,  a  typical  calculation 
with  10®  range  steps  involving  the  averaging  of  50  individual  runs  with  different  realizations  of  the 
rough  bottom  use  about  100  seconds  of  CRAY  Y-MP8  processor  time,  without  parallizing  the  code. 
This  same  calculation  can  be  done  on  a  dedicated  IBM  RS-6000/365  in  less  than  15  minutes,  and 
depending  on  the  usuage  of  the  large  scale  computer,  the  total  run  time  can  be  shorter  on  the  RS- 
6000.  A  single  3-D  run,  with  10®  range  steps  at  75  Hz  can  be  performed  in  about  7  hours  on  the 
latter  and  in  about  25  minutes  on  the  CRAY,  with  a  AO  of  20'’  without  parallelization,  which  can 
reduce  the  computation  time  by  about  an  order  of  magnitude.  The  computations  can  thus  be  seen 
to  be  practical,  but  only  on  a  supercomputer  and  only  by  parallelizing  the  code. 

3.4  Plan  of  Research 

The  development  of  the  3-D  scattering  capability  proceded  as  follows: 
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1.  Development  of  the  3-D  PE  code.  This  consisted  of  modifying  the  existing  finite-element 
code  to  output  the  complex  pressure  field  instead  of  transmission  loss.  All  of  this  development,  and 
testing,  took  place  on  a  dedicated  RS-6000,  a  RISC  cluster,  and  on  a  Convex  machine  at  Tulane 
University. 

2)  The  development  of  an  ensemble-averaging  code  which  calls  the  3-D  program,  as  in  the  2-D 

case. 

3)  A  choice  of  2-D  roughness  power  spectra  was  made,  based  on  the  extensive  literature  on 
seafloor  roughness,  [33], [34].  The  largely  geophysical  issues  associated  with  this  question  were  dealt 
with  only  in  passing,  since  these  are  not  germane  to  the  physics  of  the  problem  and  its  computational 
aspects. 

4)  Initial  computations  were  carried  out  for  a  variety  of  acoustic  environments  and  bottom 
roughness  power  spectra  as  a  test  of  the  method  and  its  implementation. 

3.5  Conclusion 

The  program  described  here  consisted  of  the  initial  stages  of  the  development,  implementation, 
and  validation  of  a  three-dimensional  finite-element  parabolic  equation  propagation  code  which 
will  have  the  capability  of  scattering  acoustic  energy  from  a  statistically  rough  bottom  by  using  an 
ensemble-average  method.  This  will  offer  the  possibility  of  accurate  modeling  in  situations  where  the 
2-D  approximation  is  inadequate.  Full  exploration  of  these  issues  will  require  significant  amounts 
of  supercomputer  time.  Such  resources  were  available  at  the  start  of  this  research,  but  became 
unavailable  in  its  later  stages.  Parallelization  is  essential,  as  discussed  above,  and  is  being  explored. 
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